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ABSTRACT 

This paper introduces a new two-parameter family of dwarf spheroidal (dSph) galaxy 
models. The mass distribution has a Plummer profile and falls like R^^ in projection 
in agreement with the star-count data. The first free parameter controls the velocity 
anisotropy, the second controls the dark matter content. The dark matter distribution 
can be varied from one extreme of mass-follows-light through a near-isothermal halo 
with flat rotation curve to the other extreme of an extended dark halo with harmonic 
core. This family of models is explored analytically in some detail - the distribution 
functions, the intrinsic moments and the projected moments are all calculated. 

For the nearby Galactic dSphs, samples of hundreds of discrete radial velocities 
are becoming available. A technique is developed to extract the anisotropy and dark 
matter content from such data sets by maximising the likelihood function of the sample 
of radial velocities. This is constructed from the distribution function and corrected 
for observational errors and the effects of binaries. Tests on simulated data sets show 
that samples of ~ 1000 discrete radial velocities are ample to break the degeneracy 
between mass and anisotropy in the nearby dSphs. Interesting constraints can already 
be placed on the distribution of the dark matter with samples of ~ 160 radial velocities 
(the size of the present-day data set for Draco). 

The Space Interferometry Mission or SIM allows very accurate differential astrom- 
etry at faint magnitudes. This can be used to measure the internal proper motions of 
stars in the nearby Galactic dSphs. Our simulations show that ~ 100 proper motions 
are sufficient to demolish completely the mass-anisotropy degeneracy. The target stars 
in Draco are at magnitudes of y 19 — 20 and the required proper motion accuracy 
is 3 — 6/iasyr~^. The measurement of the proper motions of a sample of ~ 100 stars 
uncontaminated with binaries will take about 400 hours of SIM time, or under 2% of 
the mission lifetime. 

Key words: galaxies; individual: Draco, Sculptor - galaxies: kinematics and dynamics 
- Local Group - dark matter - celestial mechanics, stellar dynamics 



1 INTRODUCTION 

The dark matter content of low-luminosity dwarf galaxies, as 
inferred from analyses of their internal stellar and gas kine- 
matics, makes them the most dark matter dominated of all 
galaxies (Mateo 1998, Carignan & Beaulieu 1989). Of these 
small galaxies, the low-luminosity gas-free dwarf spheroidals 
(dSph) are the most extreme. Available stellar kinematic 
studies provide strong evidence for the presence of domi- 
nant dark matter (e.g., Aaronson 1983, Mateo 1998), con- 
firming speculations based on estimates of the dSph's tidal 
radii (Faber & Lin 1983). Although there have been alter- 
native suggestions as to the origin of the large mass-to-light 
ratios of the dSphs (e.g., Kuhn & Miller 1989, Kroupa 1997), 
none of these have carried much conviction - see, for exam- 
ple, the objections raised by Sellwood & Pryor (1997) and 
Mateo (1997). The spatial distribution of the dark matter 
in the dSphs is very poorly known. Three obvious possibili- 
ties suggest themselves. First, the dark matter may shadow 



the stars and so the mass distribution may follow the light. 
Second, the dark matter may be distributed in a halo which 
generates a flat rotation curve, as is the case for galaxies like 
the Milky Way or M31. Third, the scale length of the dark 
matter may be larger than the luminous matter and so the 
dSph may lie in the harmonic core of an extended dark mat- 
ter halo. One promising way to distinguish between these is 
by dynamical modelling. 

Amongst the nearest dSphs are Draco and Sculptor, at 
heliocentric distances of 82 and 79 kpc, respectively. Both 
are very attractive candidates for the study of dark matter 
through dynamical modelling being relatively simple sys- 
tems with evidence for substantial dark matter content. The 
Draco dSph has an inferred central mass-to-light ratio of 
~ 60 in solar V-band units, while Sculptor has a less ex- 
treme value of ~ 10 (Mateo 1998). Any robust dynamical 
analysis is eased if the potential is approximately steady- 
state, and if the tracer stellar distribution is in equilibrium 
and well-mixed. The internal crossing times of the dSphs are 
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typically only 

tcross ~ R/(J ~ 2 X 10'^(i?/200 pc)(10 kms^Vc^) yi'- (1) 

The Draco dSph is dominated by an intermediate-age to old 
stellar population, but with little star formation for the last 
~ 5 Gyr (Hernandez, Gilmore & Valls-Gabaud 2000). Deep 
HST imaging of a small field in the outer regions of the 
Sculptor dSph has revealed a stellar population that is old 
and metal-poor, similar to halo globular clusters, with little 
spread in age (Monkiewicz et al. 1999) or metallicity (Ma- 
teo 1998). Thus, star formation indeed ceased many cross- 
ing times ago and both these systems should be well-mixed. 
Radial velocity surveys are available for both Draco (Har- 
greaves et al. 1996; Armandrofi' et al. 1995) and Sculptor 
(Queloz, Dubath & Pasquini 1995). Sculptor is of particular 
current interest due to the detection of significant amounts 
of HI gas projected within its tidal radius, at a consistent 
velocity to be truly associated with the dSph (Carignan et 
al. 1998). 

The aim of this paper is to provide new models for 
the dSphs and new techniques for probing the dark mat- 
ter distribution. Section 2 describes our models and their 
intrinsic properties, while Section 3 presents the observable 
properties, including the distributions of radial velocities. In 
Section 4, Monte Carlo simulations are used to assess how 
radial velocity surveys and proper motions inferred from as- 
trometric satellites discriminate between different dark mat- 
ter distributions. There are already samples of over a hun- 
dred radial velocities available for Draco, and this will rise 
to several hundreds in the next few years. The Space I n- 
terferometry Mission (SIM, see 'http://sim.jpl.nasa.gov/") 
has the capabilities to measure the internal proper motions 
of stars in Draco. We assess the likely impact of the new 
data sets and devise strategies for exploiting them. Finally, 
a companion paper in this issue of Monthly Notices presents 
an application of the models and algorithm to a newly ac- 
quired data set for Draco (Kleyna et al. 2001, henceforth 
Paper II) 



2 DWARF SPHEROIDAL MODELS 
2.1 Potential and Density 

In the past, dSphs have often been fitted using isotropic, 
spherical, single component King (1962) models (e.g., Hodge 
1966). These are flexible and convenient, but they do come 
at the cost of very strong assumptions - namely that mass 
follows light and that all the stars have a Maxwellian dis- 
tribution of velocities out to the tidal radius. Anisotropic 
King models relax the latter assumption, but still assume 
that the distributions of mass and light are identical. Rep- 
resentation of the dark matter and stars as two compo- 
nents of a multi-mass King model is not physically appro- 
priate as such models assume energy equipartition among 
the different mass classes, which is certainly not the case in 
a coUisionless system such as a dSph. Noting this caveat, 
Pryor and Kormendy (1990) applied such a model to data 
on the light distributions and central velocity dispersions 
of the Draco and Ursa Minor dSphs and found that mod- 
els in which the luminous and dark matter had similar dis- 
tributions were favoured over models with more extended 
dark matter distributions. However, the assumed coupling 



between the dark and luminous matter means that changing 
the velocity anisotropy of these models affects the distribu- 
tion of the luminous matter. In this paper, instead, we build 
a family of fully consistent distribution functions for the 
stars in a dSph, where we assume that the stars are tracer 
particles, moving in the underlying dark matter potential. 
This allows us to probe the mass-anisotropy degeneracy dis- 
cussed above. 

Plummer's (1911) model was originally developed to fit 
the light distribution of the globular clusters, but a much 
better application is to fit the light distribution of the dSphs 
(Lake 1990). For example. Figure 1 shows the best fitting 
Plummer profile as compared to the background-subtracted 
star count data of two Galactic dSphs, Draco and Sculp- 
tor, as given in Irwin & Hatzidimitriou (1995). With the 
exception of the Sagittarius and Ursa Minor dSphs, all the 
remaining seven Galactic dSphs are roughly spherical, with 
ellipticities lying between 0.13 (Leo II) and 0.35 (Sextans). 
So, the assumption of spherical symmetry is reasonable. We 
envisage our models as being particularly useful for Draco 
and Sculptor, because they are nearby and a wealth of kine- 
matical data is either already available or will become so 
over the next few years. 

Accordingly, let us take the luminosity density of the 
dSph as 

Po 



p{r) = 



[l + (r/ro) 



215/2' 



(2) 



where po is determined by the total observed luminosity. 
The surface brightness of the dSph is 

nR) = \ , f°^°,,2 > (3) 
3 [l + {R/T,f] 

where R is the projected radius. The physical meaning of ro 
is that it is the radius of the cylinder that contains half the 
light; henceforth ro is set to unity. 

Our aims are to assess the severity of the degeneracy 
between velocity anisotropy and mass in the dSphs and to 
investigate what radial velocity surveys may teach us. We 
need a fiexible family of dSph models with differing dark 
matter distributions and velocity anisotropies, but fitting 
the same star count profiles on the sky. So, we assume that 
the potential of the system has the form 
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[ -^log [l-f r^] ffa = 0, 
where —2 < Q < 1. Setting i/^o = v^/a, the circular velocity 
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In other words, the circular velocity curve falls off asymptot- 
ically like r'"^^ . This dark matter potential nicely spans the 
range of dark matter density distributions which we wish to 
probe: a = 1 corresponds to a mass-foUows-light Plummer 
potential and a Keplerian fall-off at large radii; a = yields 
an asymptotically flat rotation curve; and a = —2 gives a 
harmonic oscillator potential corresponding to the central 
regions of an extended dark-matter halo. Note that, as the 
parameter a decreases, the dSph becomes more and more 
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Table 1. The numerical constants in the DFs. These are pure numbers fixed once and 
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dark matter dominated. 

Figure 2 shows the variation of the central and average 
mass-to-light ratio as a function of the model parameters a 
(which controls the amount of dark matter present) and v 
(which measures the anisotropy of the velocity distribution 
- see section 2.3). It is important to bear in mind that these 
mass-to-light ratios correspond to models with the same star 
count density and the same central velocity dispersion. The 
figure provides a telling indication of the severity of the de- 
generacy between anisotropy and mass, since the mass-to- 
light ratio varies by an order of magnitude as we scan the 
models. It is this degeneracy we wish to break with kine- 
matic measurements. 



2.2 Distribution Functions 

We now build the distributions of velocities that supports 
the dSph stellar density (2) in the dark matter potential (4). 
As is well known, once the velocity distribution is permitted 
to be anisotropic, there are many possible ways of building 
a given density from stellar orbits (see Binney & Tremaine 
1987, chap. 4). The even part of the distribution function 
(DF) is determined by the stellar density law, the odd part 
by the stellar streaming or rotation law. There is no evidence 
for rotation in the Draco or Sculptor dSphs, so we calculate 
only the even part of the DF. 



2.2.1 Isotropic DFs 

According to Jeans' (1915) theorem, the DF of the stars in 
a potential is a function entirely of the isolating integrals 
of motion. Isotropic models have DFs that depend only on 

the binding energy E. The isotropic DFs are found from 
Eddington's (1916) inversion formula as: 



FiE) = 



Co,o exp (^5E/vo) if a = 0, 



(6) 



where the numerical constants Ca.o are given in Table 1. 
These are very simple results compare, for example, 
the much more complicated isotropic DFs of the spheri- 
cal isochrone (written out in Binney & Tremaine 1987) or 
the Hernquist (1990) sphere. The only other spherical mod- 
els with comparably simple DFs are the power-law spheres 
(Evans 1993, 1994). 



2.2.2 Anisotropic DFs: a > 

Anisotropic DFs of spherical models depend on both the en- 
ergy E and the norm of the angular momentum L. The DFs 
contain an additional parameter 7, which we shall see shortly 
controls the anisotropy. We begin by writing the Plummer 
density (2) in terms of the spherical polar radius r and the 
potential tp{r). There are many ways to do this, and each 
one corresponds to a different anisotropic DF. We choose 
the density partition 



Vvo/ 



5/ OL — ^/oL 



(1 + r^ 



-7/2 



(7) 



For the case of a falling rotation curve (a > 0), we can 
exploit the results in Appendix A to obtain the DF as 

F{E,L^) =Cot,^\E\^^-''^/"-'^'^H{E,L^). 

When < 2\E\, then 



(8) 
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(10) 



Here, 2F1 is the hypergeometric function of Gauss, while the 
numerical constants Ca,-y are given in Table 1. This result 
generalises an earlier calculation by Dejonghe (1987), which 
is restricted to the case where mass follows light {a = 1). 
Note that (10) is not the analytic continuation of (9) beyond 
the unit circle ~ 1- 

It is worth remarking that when 7 — — 2n, the hyperge- 
ometric function reduces to a polynomial of order n in j^, 
and then the DFs (8) are entirely elementary. For example, 
when 7 = —2, we obtain 



F{E,L^ 



\E 
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(11) 



2.2.3 Anisotropic DFs: a <0 

For the case of a rising rotation curve (a < 0), we use the 
density partition 



/ _W) \ 5/a-'r/a 



-7/2 



(12) 
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Figure 1. This shows the projected density of the best fitting Plummer model for Draco (upper panel) and Sculptor (lower panel). The 
data are the background-subtracted star count major axis profiles given by Irwin & Hatzidimitriou (1995). (ro = 9.71 arcmin for Draco 
and 13.64 arcmin for Sculptor). 



We develop the necessary formula for the anisotropic DF 
corresponding to this partition in Appendix A. The result is 

F{E,L^) = Cc,.,\Ef^^'°'-'"^G{E,L^). (13) 

When < 2\E\, the function G{E, L^) takes the form; 

G(£,L^)=2Fi(2, 1 + 2^,1;!^), (14) 



whereas when L > 2|i5|, we have; 
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(15) 
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This time, (15) is the analytic continuation of (14) beyond 
the unit circle (e.g., eq [15.3.7] of Abramowitz & Stegun 
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Figure 2. Variation of mass-to-light ratio (in units of Mq/Lq) as a function of a for three different u values. The parameter u measures 
the anisotropy of the velocity distribution. Models with u = are everywhere isotropic, while models with positive (negative) v become 
increasingly radial (tangential) at large radii (see section 2.3 for discussion). The top panel shows the total M/L within 4ro (3ro f^i Draco 
King tidal radius; Irwin & Hatzidimitriou 1995). The bottom panel shows the central M/L. 



1970). 

Again, it is worth remarking that when 7 = —2n, the 
expression reduces to a polynomial of order n in and 
then the DFs (13) are entirely elementary. For example, 
when 7 = — 2, we obtain 



F{E,L'') = Cc,-2\E\''/''-^/'' 
which is the same as above (11) 



1 _ (3 _ r^J^ 



(16) 



The DFs (18) are then entirely composed of elementary func- 
tions, namely the product of a polynomial of the angular 
momentum and the exponential of the energy. It is worth 
writing out the lowest member to emphasise its simplicity. 
When 7 = —2, we have: 



F{E,L^ 



i2nvl) 



1 + 



71/ 
2^ 



exp 



(20) 



2.2.4 Anisotropic DFs: a = 

The case of a flat rotation curve {a = 0) has a different form 
again. We write the Plummer density as 

p = po exp ((5 - 7)V'/«o) (1 + r')-^/'. (17) 

The DF corresponding to this density partition is derived in 
Appendix A. We find: 

F{E, L') = Co,, exp ((5 - ^)E/v'o) 1, ^^^) , (18) 

where $ denotes the degenerate hypergeomotric function 
and the constant Co, 7 is given in Table 1. When 7 = — 2n, 
the function $ reduces to the Laguerre polynomial L„ (see 
eq. [8.972.1] of Gradshtcyn & Ryzhik 1978): 



-{2n + 5)L^ 
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exp^(2n-F5)£;/wo)- 



(19) 



2.2.5 Practical Evaluation of the DFs 

The DFs (8), (13) and (18) all depend on the hypergeomet- 
ric function in one way or another. Numerical evaluation of 
the hypergeometric function can be difficult. In fact, none 
of the algorithms for evaluating the hypergeometric function 
2Fi{a,b, c; z) presented in the standard reference books on 
numerical methods (e.g.. Press et al. 1992) is applicable for 
all values of the parameters a, b and c. This is because very 
different strategies are required depending on the magnitude 
and signs of the parameters. Our computational algorithm 
is presented in Appendix B. Note that in our applications to 
calculations of the likelihoods in Section 4, it is important 
that the DFs are calculated extremely accurately. The obvi- 
ous way of checking our numerical algorithm is to investigate 
whether the integration of the DF over velocity space yields 
the Plummer density. Typically, we find that the density is 
recovered to better than one part in lO". 
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2.3 Intrinsic Moments 



3 OBSERVABLE PROPERTIES 



Although the DFs (8), (13) and (18) have differing forms ac- 
cording to whether the rotation curve is falhng (a > 0), flat 
(a = 0) or rising (a < 0), nonetheless all physical quantities 
like the moments vary smoothly and continuously with the 
rotation curve index a. 

The intrinsic velocity second moments are 
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(22) 



The second moments converge and are positive definite pro- 
vided 7 < min(a-|-5, 2). 

In terms of Binney's (1981) anisotropy parameter a, the 
radial and tangential velocity dispersion {v^) and {vg) vary 
as 

„2 
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(23) 



If 7 = 0, the velocity dispersions are isotropic. Irrespective 
of 7, the central regions of the model are always isotropic. 
At large radii, the anisotropy becomes constant (c.f , Henon 
1973; Wilkinson & Evans 1999). If 7 < 0, then the dispersion 
tensor becomes tangential with increasing radius; if 7 > 0, 
it becomes radial. The limit 7 —00 is the circular orbit 
model, while the limit 7 — > 2 is built from radial orbits alone 
in the outer parts. 

It is also helpful to define the quantity v which is related 
to 7 by 



logi 



■ 7 



This quantity runs from —00 (tangential velocity distribu- 
tion at large radii) to -|-oo (radial velocity distribution at 
large radii). The velocity distribution is everywhere isotropic 
for u — 0. This definition is helpful because the ranges of 
tangential and radial anisotropy are symmetric about v = Q. 

The intrinsic fourth moments are useful for diagnosing 
deviations from pure Gaussianity and so we briefly list the 
results here: 
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(25) 



The properties of dSphs when viewed in projection are of 
particular interest to us, since these correspond to observ- 
able quantities. In this section, we present the distributions 
of the line of sight velocities and the proper motions. All the 
moments of these distributions are analytic. 



3.1 Line of Sight Moments 

To derive the line of sight velocity moments involves per- 
forming a triple integration along the line of sight and over 
the two velocity components perpendiculax to the line of 
sight. For our dSph models, these integrals are all analytic. 
The line of sight second moment is given by 



2 



{1 + R 
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7(4+q) 



2(5 + a) l + i?2 



where the central velocity dispersion is 
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3V^t;gr(2-t-a/2) 
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(27) 



(28) 



When a = 1, this reduces to the result found by Dejonghe 
(1987) for his self-consistent Plummer models. The left pan- 
els of Figure 3 show the variation of the line of sight sec- 
ond moment with projected radius R for a = —1,0, 1 and 
7 = 2,0,-10 (or, equivalently, v = (X),0, — 0.8). For a > 
and 7 < — Q!(5+a)/(4r|-a), the curves (7p{R) show a maximum 
at 



^2 ^ 2(Q(5 + a) + 7(4 + a) 
a(7(4+Q) - 2{5+a)) ' 



(29) 



As Figure 3 indicates, the models have an unusual and at- 
tractive feature. For a given a, all the curves pass through 

the same point 
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irrespective of 7. This radius is the changeover radius. Ra- 
dially anisotropic dSphs (7, v > 0) have less projected ve- 
locity dispersion ap at large radii than isotropic dSphs, and 
more ap at small radii. This follows from simple geometrical 
considerations. Radial orbits contribute significantly to the 
line of sight motion in the inner parts, but much less so in 
the outer parts. Conversely, tangentially anisotropic dSphs 
(7, z/ < 0) have more projected velocity dispersion crp at 
large radii than isotropic dSplis, and less ap at small radii. 
The changeover radius is the radius at which the transition 
from the central to the asymptotic properties occurs. 

The line of sight fourth moment is useful for detecting 
deviations from Gaussianity. It can be calculated from the 
intrinsic fourth moments as: 



{ve) = {v<i>) = 
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All the remaining components of this fourth rank tensor van- 
ish. The fourth moments converge and are positive definite 
provided 7 < min(2Q-|-5, 2). 
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with the central value given by 
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Figure 3. This shows the Vciriation of the line of sight second moments (left panels) and kurtosis (right panels) with radius R on the 
plane of the sky for models with different dark matter content (a = 1,0, —1 as labelled). The full curves refer to the isotropic model 
{ly = 0), the dashed curve to a radially anisotropic model {u = oo), the dotted curves to a tangentially anisotropic model (f = —0.8). 



In fact, all the line of sight moments are analytic. The nth 
moment CTp" is 

^p"(^) = (TTil)w^ (-^' i' 1' w) '(33) 

where 3F2 is the generalised hypergeometric function (which 
always reduces to a finite nth order polynomial). The central 
value is 

2„ _ ^n-2 3V.gr(n+l/2)r(5/a-7/a+l)r(2+na/2) 

r(5/a-7/a+n+l)r(5/2+na/2) ' ^ ' 

Figure 3 shows the variation of the kurtosis k — 
(jp/ia^Y with projected radius R. The kurtosis measures 



the extent to which the distribution is peaked. A Gaussian 

distribution has a kurtosis of 3. The larger the kurtosis, the 
broader is the distribution. The isotropic models (7, v = Q) 
always have a constant kurtosis given by 



4 r(Q+2)r(Q/2 + 5/2)r(Q/2 + 7/2) 
v/¥ r(Q + 7/2)r2 (a/2 + 2) 



(35) 



In fact, the a = 7, = model has a constant kurtosis of 3, 
indicating that the line profiles of this model are Gaussian. 
Figure 3 illustrates the tendency for increasing kurtosis as 
the dSph becomes more dark matter dominated. This stems 
predominantly from the larger tails in the intrinsic velocity 
distributions. As a diminishes, there are more and more high 
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Figure 4. Normalised line profiles or distributions of line of sight velocities fios are shown for models with different dark matter content 

(a = 1,0,-1 as labelled) and different anisotropy {y = —0.6,0,0.6). The line profile is calculated at i? = 0, the center of the dSph in 
projection. For comparison, the distributions of proper motions {vr and v^) on the plane of the sky are also shown - these profiles are 
identical for R = Q and appear superposed as a dot-dashed line. 



velocity stars in the tails. 
3.2 The Line Profiles 

When comparing models to the discrete radial velocities, 
the main quantity of interest is the line profile. Tiiis is the 
probability distribution of the line of sight velocities at a 
given projected radius. Mathematically, the unnormalised 
line profile L(wios, R) is the integral of the DF along the line 
of sight and over the tangential components of velocity. Let 
us use {R, ip) as polar coordinates on the plane of the sky 
and z as the coordinate along the line of sight. 
When a < 0, the line profile is: 



dz / dvR / dv^F{E,L^), 

-OO J —OC J —OC 



(36) 



whereas when a > 0, the line profile is: 

L{v\as,R) 



dz I dvR 
X / dv^F{E,L^). 

(2-^^ _„2)l/2 
^ los 



(37) 



If Q < 0, then the integration over the line of sight extends 
indefinitely. If a > 0, the stars with velocities vios are seen at 
projected position R only if vf^^ < 2ip{R,z). This provides 
upper and lower limits z± to the line of sight integral. Often 
it is useful to divide by the surface brightness S(i?) and con- 
sider normalised line profiles £{v\ob, R) = L{vios,R)/^{R)- 

For the isotropic models, the line profiles cam be reduced 
easily to single quadratures. When a ^ 0, the unnormalised 
line profiles have the form 



L{vioa, R) 



2 KaPO 



f 

J R 



b/a-l/2 



rdr 



(38) 



where 

r r(5/a+l) 
T{5/a+l/2) 



r(l/2-5/a) 
r(-5/a) 



(39) 



if a > 0. 



If and only if the rotation curve is flat (a — 0) and the model 
is isotropic (7,1/ = 0), then the normalised line profile is a 
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Figure 5. As Figure 4, but for locations at projected position R = 4. 



"1 ' r 



J I L 



"1 I r 



pure Gaussian 



{vios,R) 



(4)'"-(-ir) 



(40) 



We emphasise that the Hne profiles vary smoothly and con- 
tinuously through the seemingly singular case of a = 0. This 
may be verified by taking the limit a ^ of (38). 

Although presently unmeasurable with ground-based 
telescopes, the internal proper motions of stars in the nearby 
dSphs are accessible to future space missions (like the Space 
Interferometry Mission or SIM). In Draco and Sculptor, the 
stars of interest are at magnitudes V ~ 19 ^ 20 and tlic re- 
quired proper motion accuracy corresponding to 1 — 2 kms~^ 
in velocity is 3 — 6/Ltasyr~^, which is just within the ca- 
pabilities of SIM. Hence, it is also useful to calculate the 
one dimensional distributions of proper motions L{vr,R) 
or L{vtp,R) by integrating along the line of sight and over 
the other two transverse velocities. 

Figures 4 and 5 present some examples of line profiles 
and proper motion distributions obtained for a variety of val- 
ues of R, a and 7. These are calculated by direct integration 
over the DFs. For the isotropic models — 0), the line 

profile and the proper motion distributions clearly coincide. 
Tangentially anisotropic clusters tend to show bimodal line 
profiles, whereas radially anisotropic clusters tend to have 
narrower, peaked line profiles. Hence, the trend as we move 



vertically downwards in the panels is towards broader and 
more flat-topped profiles. The larger the dark matter con- 
tent, the greater and more distended the wings of the line 
profiles. So, the trend as we move left to right in the panels 
is towards less extended line profiles. 

Note that, as we move outwards from the center of the 
dSph, the differences between the line profiles become more 
pronounced. Observationally speaking, it is easier to obtain 
radial velocities for stars in the centre, but it is stars in the 
outer parts of the dSph which are most useful for breaking 
the degeneracy between mass and anisotropy. 



4 METHODOLOGY 

Previous studies of dSph dynamics have been impeded by 
the degeneracy between anisotropy and mass. An increase 
in the line of sight velocity dispersion at large radii may 
by due to either (1) the presence of large amounts of mass 
at large radii, or (2) tangential anisotropy in the velocity 
distribution. The large observed central radial velocity dis- 
persion is compatible with either a massive halo, and a low 
central density; or no halo, and a large central density. In 
this section, we test our ability to discern between these pos- 
sibilities using (i) large radial velocity surveys and (ii) radial 
velocities and proper motions provided by SIM. 
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Figure 6. Distribution of binary velocities as inferred from the binary distribution in the solar neighbourhood. Dashed curve - 
distribution. Solid curve - evolved distribution for giant star binaries. See text for details. 
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4.1 The Likelihood Function 

The problem of comparing discrete radial velocity measure- 
ments to galactic models has already received much atten- 
tion (e.g., Little & Trcmaine 1987, Kochanok 1996, Wilkin- 
son & Evans 1999, Evans & Wilkinson 2000). Suppose wc are 
given radial velocities Viosi (corrected for the galaxy mean 
velocity) of N stars in a dSph at projected positions Ri. At 
each point, we can compute the probability of observing the 
radial velocity data set 



;a,7)- 



(41) 



So, we scan a grid of a, 7, and at each point compute the 
probability of observing the entire input data set. Using 
Bayes' theorem and assuming uniform prior probabilities 
in the model parameters, then the most likely values of a 
and 7 are given by maximising (41). Confidence regions axe 
obtained by applying two-dimensional statistics to the 
likelihood (or the logarithm of (41)). 

This procedure requires repeated evaluation of the line 
profiles, and direct integration over the DF is too slow and 
expensive to provide a competitive algorithm. Accordingly, 
we only use brute force integration to provide a look-up 
table of the line profiles in a grid in the four-dimensional 
(q,7, 7?,w) space. The grid spacing in a, R and v is linear, 
whereas the spacing in 7 is uniform in log(2 — 7). We use cu- 
bic splines to interpolate in the logarithm of the line profile 
between the grid points. Once the look-up table has been 
built, this provides an extremely fast and accurate way of 
calculating the line profiles. Typically, the error in the line 
profile (as inferred by integration over the line of sight ve- 
locity to recover the surface density) is better than one part 
in lO"*. This is still sufficient to mislead our maximum like- 
lihood algorithm because the probability (41) is formed by 



the multiplication of A'^ line profiles. To take account of this 
small error, we re-normalise each line profile to unity after 
interpolation. 

Before applying the likelihood algorithm, there are two 
further corrections that must be applied. First, the line pro- 
file is convolved with a Gaussian of width 2kms~^ to allow 
for observational errors. This is a typical error for data ob- 
tained on the 4m class telescopes like the William Herschel 
Telescope with multi-object spectrographs. Second, the line 
profile is adjusted for the efi^ect of binaries. This also involves 
convolution, this time with the binary correction function 
h{v). By Monte Carlo sampling binary orbits drawn from 
the binary distribution in the solar neighbourhood, taking 
account of tidal circularisation of the orbits as described in 
Paper II, we deduce the distribution of velocities Pb{v) in- 
duced by the binary motion. This is shown in Figure 6. The 
binary correction function is then given by 

h{v) = fP^{v) + (1 - /)5(«), (42) 

where the constant / is the binary fraction. 

4.2 Radial Velocity Surveys 

Synthetic data sets are created by choosing phase 
space coordinates {Ri, Zi,viii,Vip^,viosi}i=\...N drawn from 
the DF, and then discarding all but the projected positions 
and the line of sight velocities. The velocities are contam- 
inated by measurement noise of amplitude 2kms~^, and a 
fraction / of the stars are also assigned a binary velocity 
which is added to their line of sight velocity. To reduce the 
computation time, rather than randomly selecting the pro- 
jected radius for each star, ten values of the projected radius 
-Ri between 0.0 and 2.5 are used. Each data set is then anal- 
ysed and a most likely value of a and 7 (or, equivalently, v) 
is obtained. 
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Figure 7. The recovery of the dark matter parameter a and the 
anisotropy u from simulated data sets of 160 radial velocities. Ob- 
servational errors of 2 km s~^ are assumed for the radial velocities 
and a binary fraction of 40% is assumed. The main panel shows 
where each data set falls in the {a, v) plane. The side panels show 
the one-dimensional histograms. The true values of a and v are 
marked with arrows. 
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Figure 9. The recovery of the dark matter parameter a and 
the anisotropy v for simulated data sets of 160 stars with 
ground-based radial velocities, together with proper motions mea- 
sured by SIM. The true values of a and v are marked with arrows. 
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Figure 10. As Figure 9, but for simulated data sets of 1000 stars 
with radial velocities and proper motions. 



Figure 8. As Figure 7, but for simulated data sets of 1000 radial 
velocities each. 

Figure 7 shows the results from the analysis of simulated 
data sets of 160 radial velocities. This is the number of radial 
velocities that are presently available for Draco. Figure 8 
shows the results from analysis of simulated data sets of 1000 
radial velocities. This is an estimate of the maximum number 
of radial velocities that may be available in Draco with an 



observing program on a 4m class telescope. In each figure, 
the scatter plot shows the joint distribution of a and v values 
obtained from the artificial data sets while the histograms 
show the spread in the individual parameter values. 

Even with a data set of the present size, Figure 7 shows 
that it is possible to put interesting constraints on the value 
of a. The input model for Figure 7 had an ct value of —0.5 
indicating a dark matter distribution intermediate between 
that of a flat rotation curve halo and an extended harmonic 
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Figure 11. The recovery of the dark matter parameter a and the 
anisotropy v for simulated data sets of 160 stars with radial veloc- 
ities. Data produced with a binary fraction of 80% but analysed 
assuming a binary fraction of 40%. 




Figure 12. The recovery of the dark matter parameter a and 
the anisotropy v for simulated data sets of 160 stars with radial 
velocities. Data generated from a model with a Jaife sphere rep- 
resenting the dark halo of the dSph. See text for details. 



core. From the figure, the 2a confidence interval on a is 
(—1.45,0.8). Thus data sets containing 160 radial veloci- 
ties are sufficient to rule out a mass-foUows-light halo model 
(a = 1) at about the 2.5cr level. Increasing the size of the 
data set to 1000 stars allows much greater discrimination 
between models. 



4.3 Proper Motion Surveys 

Here, we consider the possible impact of next generation 
astromctric satellites like SIM and its successors. Microarc- 
second astrometry to an accuracy 3 — 6/Ltasyr~^ is sufficient 
to allow the internal proper motions of the brightest stars in 
Draco to be measured to 1— 2kms~'^. Since all stars in Draco 
are at almost the same distance, only the relative proper 
motions arc required and this significantly reduces the ob- 
serving time. A further advantage is that all the target stars 
are within a ~ 1° field of view, allowing subdivision of the 
observations into a small number of overlapping fields and 
hence reducing the instrumental systematic contribution to 
the error. A satellite like SIM is ideal for these measure- 
ments, as it allows very accurate differential measurements 
at faint magnitudes. 

To simulate this, we pick phase space coordinates 
{xi,yi,Vxi,Vy^,Vzi}i=i...N drawn from the DF, discard the 
geocentric radial position coordinate z , and attempt to re- 
cover a and v using the Bayesian likelihood. As in the rar 
dial velocity simulations described above, we contaminate 
the data with measurement noise in each component of ve- 
locity. The correct expression for the probabilities required 
to analyse these data is then the 3-dimensional convolution 
of the DF with 3 Gaussians, each of width 2kms~^. How- 
ever, as was noted earlier, integration of the DF is com- 
putationally expensive, and we therefore approximate the 
convolution by a weighted sum over the values of the DF at 
the corners of a cube of side \/2cr. This approach gives suffi- 
ciently accurate results for the present paper - in the future, 
increases in computational speed will allow us to perform the 
full 3-dimensional convolution. At present, wc consider only 
the case where no binaries are present. This is both com- 
putationally convenient, and may be justified on the basis 
that proper motion measurements will be based on multi- 
ple epoch observations which will allow the identification of 
binaries in our sample. 

Figure 9 shows the results from the analysis of simu- 
lated data sets of 160 with both radial velocities and proper 
motions measured for each star to an accuracy of 2kms^^. 
Figure 10 shows the results from analysis of similar data sets 
containing 1000 stars. The narrow spread visible in these his- 
tograms, particularly in the case of 1000 stars, emphasises 
the value of having all components of the stellar velocities 
when modelling the mass profile. The key feature of these 
results, as illustrated in Figure 9 , is that even with only 
160 proper motions it is possible to break unambiguously 
the degeneracy between mass and anisotropy. 



4.4 Robustness 

In this section, we consider the robustness of our analysis to 
errors in the estimation of the sample binary fraction and 
measurement errors. As the same algorithm is used in both 
cases, we present results only for the former case. We also 
address the question of the extent to which our analysis is 
model dependent. 

Let us begin by considering data sets in which the bi- 
nary fraction has been severely underestimated. Data sets 
of 160 radial velocities are generated as described earlier, 
assuming a binary fraction of 80%, and measurement errors 
of 2kms~^. However, when analysing the data, we perform 



Dwarf Spheroidals I: Models 13 



the convolution of the hne profiles assuming a binary frac- 
tion of only 40%. This discrepancy represents an upper limit 
to the uncertainty in the binary fraction in a typical data 
set of single-epoch radial velocities. The results of this ex- 
periment are shown in Figure 11. Comparison of this figure 
with Figure 7 shows that while the distributions of recovered 
Q and 7 values are slightly broadened, the overall effect of 
the mis-estimation of the binary fraction is not serious. The 
uncertainty caused by the true size of observational errors 
can be checked in a similar manner. Wc generate data sets 
of 160 radial velocities with a binary fraction of 40% and 
measurement errors of 4kms~^. However, we convolve the 
line profiles with an error Gaussian of width 2kms~^ rep- 
resenting an over-optimistic estimation of the measurement 
errors. We do not show a figure for this case, as the results 
are predictable. Our algorithm is not confused and the re- 
sults are clustered around the correct values, albeit with a 
greater spread. 

A more serious question is whether our maximum like- 
lihood algorithm - which assumes a parametrised fit for the 
density and potential - is robust, as any dSph undoubtedly 
deviates to a greater or lesser extent from the model. To test 
this, we generate synthetic data from a very different halo 
model and analyse them using our models. We retain the 
assumption that the light in the dwarf follows a Plummer 
profile, but we replace the halo with a Jaffe (1983) sphere. 
The halo potential is then given by 

ip{r) = Vclog (— ^) 

where Vc is the amplitude of the halo rotation curve at small 
radii, and r.j is the scale radius of the halo. For these simu- 
lations, we set rj = 1 and Vc = 45.5 kms~^, yielding a halo 
whose rotation curve falls off for r > 1 but for which total 
mass enclosed within r ~ 3 is similar to that of an a ~ —2 
(rising rotation curve) model from our family of Plummer 
halo models. Following Wilkinson & Evans (1999, Eq. (19)) 
we build a DF (with constant velocity anisotropy P) of the 
form F{E,L) = L''^^ f{E). The value of /3 is set to be the 
mean value of the velocity anisotropy between r = and 
r = 3 for a model with v = —0.2 from our Plummer family. 
We analyse the data using radial velocities alone. 

Figure 12 shows the results obtained from data sets of 
160 radial velocities. The analysis of the Jaffe model data 
sets returns Plummer model parameter values in the region 
of a = —1.5, v — —0.2. Qualitatively, this result is very en- 
couraging, as the parameter estimates all lie in the region 
of the (a, 7) plane where the dark matter is more extended 
than the light in the dSph. More quantitatively, the Jaffe 
data sets return approximately the correct value of f. The 
mass enclosed within r = 3, the region probed by the gener- 
ated data points, is also recovered to better than a factor of 
two - for the (input) Jaffe model we obtain 3.6 x 10® M0, 
while for the recovered Plummer halo we obtain 2.2 x 10* 
M0. We coiiehide that the parameter values returned by 
our method are reasonably robust, even for a small data 
set, and that they do reflect underlying physical properties 
of the data. 

5 CONCLUSIONS 

This paper has presented a flexible family of spherical mod- 
els suitable for representing dwarf spheroidal (dSph) galax;- 



ies. The models have a Plummer profile, which is an excellent 
fit to the star count data on the nearby dSphs. There are 
two free parameters. The first is the dark matter parameter 
a. When a = 1, the mass follows the light and the rotation 
curve is asymptotically Keplerian; when a = 0, the dark 
matter is distributed in a cored isothermal sphere and the 
rotation curve is fiat; when a = —1, the dSph is enclosed 
within the harmonic core of a much larger daxk matter halo. 
The second parameter is the anisotropy parameter 7. When 
7 > 0, the models are radially anisotropic at large radii; 
when 7 = 0, they are isotropic; when 7 < 0, they are tan- 
gentially anisotropic at large radii. 

These models offer considerable advantages over King 
models, which are conventionally used for fitting dSphs 
(see e.g. Binney & Tremaine 1987). The single component, 
isotropic King models often used for modelling the dSphs 
assume that the mass follows the light and have an isotropic 
stellar velocity distribution. The cost of such strong assump- 
tions is that they are unable to probe the degeneracy that 
exists between dark matter mass and anisotropy. A very at- 
tractive property of our new models is that the intrinsic and 
projected moments are simple and analytic. In particular, 
the line of sight velocity dispersion is available as a function 
of projected radius for all values of the dark matter and 
anisotropy parameters. The phase space distribution func- 
tion (DF) is more complicated - as it depends on transcen- 
dental functions - but it is feasible to evaluate numerically. 
Quadratures over the DF provide us with the distributions 
of radial velocities (the line profiles) as well as the distribu- 
tions of proper motions. All these observable quantities are 
readily available for our family of dSph models! Note that 
the line profiles of all the models are similar in the centre, 
but begin to differ at ~ 2ro, emphasising the importance of 
gathering data at large radii. 

We have used our new models to assess what largo sam- 
ples of radial velocities can tell us about dark matter in 
dSphs. All of the information on the kinematics is contained 
in the DF. So, given large samples of stars with projected 
positions and radial velocities, we calculate the likelihood of 
observing the data set as a function of the dark matter and 
anisotropy parameters. The distribution of radial velocities 
is corrected first for observational errors by convolution with 
a Gaussian and second for the contamination by binaries by 
convolution with a suitable correction function. 

Simulated data sets of 160 and 1000 radial velocities 
are drawn from the DF, contaminated with measurement 
noise and binary velocities, and fed into the likelihood algo- 
rithm. A sample size of hundreds of velocities is typical of 
the presently available data sets, whereas a sample size of a 
thousand velocities is an estimate of the number that may 
be accessible to a long-term observing program in the nearby 
dSphs like Draco. Our simulations show that a sample size 
of 160 stars is already sufficient to make interesting state- 
ments about the dark matter distribution in a dwarf galaxy. 
In particular, using data drawn from a model with a slightly 
rising rotation curve it is possible to rule out mass-follows- 
light models at about the 2.5a confidence level. Larger radial 
velocity data sets allow much tighter constraints to bo placed 
on the model parameters. Our machinery is set to work in a 
companion paper (Kleyna et al. 2001) to interpret a newly 
acquired data set for the Draco dSph. 

The addition of proper motion data completely breaks 
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the degeneracy between mass and velocity anisotropy, even 
using only a hundred or so proper motions. Obtaining such a 
data set will be feasible using the SIM satellite. One strategy 
is to construct a sample of stars that are uncontaminated by 
binaries before SIM flies, for example, by taking second and 
third epoch data on ground-based telescopes and removing 
all radial velocity variables. SIM time is very expensive and 
it is wasteful to spend it following the astrometric paths of 
binaries if our aim is to measure the proper motions of single 
stars. In Draco and Sculptor, the target stars are at mag- 
nitudes 1/ ~ 19 and the required proper motion accuracy 
corresponding to 1 — 2 kms~^ in velocity is 15 — 30 /^as over 
the mission lifetime of 5 years. The optimum strategy is to 
measure the positions of the stars twice, once at the begin- 
ning and once at the end of the mission. This will be ample 
to calculate the proper motions, provided all binaries have 
been eliminated by our ongoing radial velocity surveys. SIM 
takes ~ 1 hour to measure the position of ~ 19 star to 
~ 20 /^as (in one dimension), which is our typical required 
accuracy. Hence, for each star, we require 4 hours of SIM 
time (for two-dimensional measurements at the beginning 
and end of the mission). Our simulations have shown that 
samples of ~ 100 proper motions are ample to break the 
mass-anisotropy degeneracy. This will take about 400 hours 
of SIM time, or roughly 1.7% of the mission lifetime. In other 
words, this is an extremely competitive use of SIM time. 

Given the modest amount of SIM time required, our 
collaboration is considering an ambitious program that ac- 
quires the proper motions of roughly a hundred stars in each 
of Draco, Ursa Minor and Sextans. This program will still 
only consume 5% of the mission lifetime. Such a data set, to- 
gether with the sophisticated modelling techniques we have 
introduced in this paper, will make it possible to map the 
dark matter distribution in these three dSph galaxies with 
unprecedented accuracy. This important test of the nature 
of dark matter on such small scales is uniquely possible with 
SIM. 
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APPENDIX A: DERIVATION OF THE 
DISTRIBUTION FUNCTIONS 

First, let us summarise Dejonghe's (1986) results, which ap- 
ply to the models with falling rotation curves (a > 0). Let 
us start with the well-known result (Fricke 1951, Dejonghe 
1987, Evans 1994) that the density p' 

I IP -a 



corresponds to the DF 

r(p+i) 



(43) 



(44) 



r(l-/?/2)(p-l/2-|-/3/2) (27r)3/22-'3/2 ' 

These are sometimes called Fricke components. From this 
simple result, we can build up the DF F corresponding to 
the density 

^2a 



(l + r2)-+''' 



(45) 



by a continuous superposition of the Fricke components. The 
result is 

V{p+l)E^-^''^ 



F{E,L') = 



where 



(27r)»/2r(a -I- 6) 



H{a,b,p-\,1-^), (46) 



H{a,b,c,d-x) = — / ^ , 'tX 'ds. (47) 



2TTi J„ T{c + s)T{d- sY 
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Here, C is a contour in the complex s plane such that all the 
poles ~-a~n are on the left and all the poles b + n are on the 
right (where n is an integer). Dejonghe (1986) shows how to 
evaluate this integral. When a; < 1, the contour is completed 
to the left by adding a large semi-circle at infinity. Only the 
poles at s = —a — n are enclosed, giving 

n{a, b, c, d; x) = /"^"'^J']^ ,. 2^1(0+6, l+a-c, a+d; x). 
i (c— a)i (a+d) 

When a:: > 1, the contour is completed to the right and only 
the poles at s = 6 + n contribute to give 

Hia, b, c, d; x) = f^I^)f^:p^2Fi(a+6, 1+b-d, b+c; i). 

Second, let us extend this result to the models with 
rising rotation curves (a < 0). In this case, the Fricke com- 
ponents are derived by Evans (1994) as 



_ r(p - 13/2 + 3/2) 



r(l - /3/2)r(p) (-£;)p-/3/2+3/2 (27^)3/22-/3/2 ' 

Again, we seek the DF F corresponding to the density 



P = 



This time the result is 



F{E,L') = 



1 



(27r)3/2r(p)r(a + 6) 
1 



(_£;)p+3/2 



6;(a,6,p-F|,l;|^), 



where 

Q{a, h, c, d; x) 



J_ f r{a+ 
2m J ^ 



s)T{b-s)V{c-s) 

r(d-s) 



(48) 

(49) 

(50) 
(51) 

(-a;)"''ds.(52) 



Here, C is a contour in the complex s plane such that all the 
poles —a — n are on the left and all the poles b + n and c + n 
are on the right (where n is an integer). When \x\ < 1, we 
obtain 

^/ L J X T{a+b)V{a+c)(-xY ^ , , 

g{a, b, c, d; x) = Y^^^i+d) {a+b,a+c, a+d; x). 

When \x\ > 1, the contour is completed to the right and 
there are two infinite sequences of poles at 6 -|- n and c + n. 
We therefore obtain a sum of two hypergeometric functions, 
namely 

g{a, 6, c, d; x) =^i^Z^I^±^,Fi{a+b, 1+6-d, l+&-c; i) 
r(b-c)T(a+c) ^ , , , , , 1 , 

+ (-.)cr(rf-c)^ ^^('^+^' 

The third and final case is that corresponding to a flat 
rotation curve (a = 0). The elementary Fricke components 
become 



p =r '^^^exp{p^/vo), 
and 

„3/2-/3/2 

/' 



(27r)3/22-'3/2r(l - /3/2)i;^ 



3-/3/2 



(53) 

L-'^/^exp{pE/vo). (54) 



This means that the density 

»,2a 

: e:xp{pip/vo), 



1^ (l-|-r2)»t* 
corresponds to the DF 

„3/2+a 



(55) 



F = 



P 



■ eMpE/voMa}i,aH,-^),i56) 



(27r)3/22ar(a-hl)wg"+3 
where $ is the degenerate hypergeometric function. 



APPENDIX B: NUMERICAL COMPUTATION 
OF HYPERGEOMETRIC FUNCTIONS 

The hypergeometric function is defined within the unit circle 
a < 1 by the hypergeometric series 



2Fi{a,b;c;z) = IH —z + 

c 1! 



a{a+l)b{b+l) 2 
c(c+ 1)2! ■ 
a{a+l) . . . ja+n-l) b{b+l) . . . jb+n-l) 
c{c+l) . . . (c-l-n-1) n! 



(57) 



As the coefficients of the series approach 1 as n — » 00, the 
series is bounded by a majorizing geometric series, and con- 
verges for \z\ < 1. Outside of the unit circle, an analytic 
continuation of the hypergeometric series is given by the 
hypergeometric differential equation 



2F1 = 



ab2Fi - [c-{a+b+l)z]2Fi 



(58) 



z{l-z) 

A second, linearly independent solution to the hypergeomet- 
ric differential equation is z^~'^ 2^^1(0— c-|-l, b—c+1; 2 — c; z) 
(provided c is not an integer). 

Numerical evaluation of the hypergeometric function 
can be awkward. In fact, none of the standard reference 
books on numerical methods (e.g.. Press et al. 1992) present 
completely general algorithms for the computation of the hy- 
pergeometric function valid for all parameter values. Luke 
(1977) has derived several useful Chcbyshev and rational 
polynomial approximations. The coefficients of the hyperge- 
ometric series (and its approximations) may reverse sign, so 
that the final value can depend on the near-cancellation of 
very large terms, rendering the normal fifteen digit floating 
point precision of a computer insufficient. The parasitic sec- 
ondary solution of the differential equation, by virtue of its 
steep z^~'' growth, renders the integration of the differential 
equation numerically impracticable for c <C 0. 

Our algorithm divides the evaluation of the hyperge- 
ometric function into two regimes. Provided c > —3, the 
hypergeometric differential equation (58) is stable. It is 
integrated via the standard ordinary differential equation 
integration methods provided by Press et al (1992). The 
initial condition for the integration cannot be chosen as 
2-Fi(a, fe, c; 2) = 1 when z = 0, as the differential equation 
(58) is then singular. Accordingly, we start from z = e and 
calculate the value of 2F1 [a, b, c; z) at this point by direct 
summation of the hypergeometric series. This is the method 
of choice for evaluation of the DF when a > or when a < 
and |LV(2S)| < 1. 

If c < — 3 (as occurs in the expressions for the DF when 
a < and \L'^/{2E)\ > 1) and z > —0.5, we compute the 
power series using standard double precision arithmetic, but 
switch to arbitrary precision arithmetic whenever the terms 
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of the series exceed 10 , ensuring that the final value will be 
accurate to nine decimal places. Arbitrary precision arith- 
metic is necessary only for a — > 0, or c <C 0. If « < —0.5 then 
we use the transformation (see eq. [9.131.1] of Gradshteyn 
& Ryzhik 1978) 

2F^{a,h-c-z) = (l-z)-"2Fi(a,c-&;c;2/(«-l)) (59) 

to cast z into z/{z — 1), ensuring that the power series con- 
verges swiftly for all z of interest. 

This paper has been produced using the Royal Astronomical 
Society /Blackwell Science 1^ macros. 



